Detection of crop lines with Hough Line Transform¶

image.png

The Hough transform for line detection is a popular technique for identifying lines in images, including drone images. The Hough transform converts image points into a parametric representation of lines, allowing the detection of lines even in images with noise or discontinuities.

*Steps*

  • Image preprocessing: Before applying the Hough transform, it is useful to carry out some image preprocessing steps to improve the quality of the results. This can include the application of filters to reduce noise, such as a smoothing filter (like the Gaussian filter) or an edge detection filter (like Canny).

  • Edge detection: The Hough transform operates on the edges of the image, so it is necessary to apply an edge detection algorithm to the preprocessed image. The Canny algorithm is widely used for this purpose, as it is effective for detecting edges and minimizing noise response.

  • Application of the Hough transform: The Hough transform is applied to the edge image to identify the parameters (angle and distance) of the lines present in the image. The standard Hough transform operates in the space of angle and distance parameters (it can be expanded to detect lines with additional parameters such as thickness).

  • Threshold and line detection: after applying the Hough transform, you must define a threshold to identify significant lines in the image. It is considered that the points in the Hough transform that surpass the threshold represent the detected lines. You can adjust the threshold to control the detection sensitivity.

  • Draw the lines detected in the original image: Finally, based on the parameters of the lines detected in the previous step, you can draw the lines that are found in the original image. This will allow you to see the identified plantation lines.

It is important to mention that the Hough transform can be sensitive to different lighting conditions, textures and other factors present in the images. Therefore, it may be necessary to adjust the parameters of the Hough transform and carry out tests in different scenarios to obtain the best results.

Furthermore, consider that detecting crop lines in drone images can entail specific challenges, such as the presence of vegetation, lighting variations and superpower lines. Therefore, it may be necessary to explore additional techniques, such as image segmentation or custom filters, to address these situations and improve the accuracy of crop hill detection.

We start by installing the Rasterio.

In [ ]:
!pip install rasterio

Requirement already satisfied: geopandas in /usr/local/lib/python3.10/dist-packages (0.13.2)
Requirement already satisfied: fiona>=1.8.19 in /usr/local/lib/python3.10/dist-packages (from geopandas) (1.9.5)
Requirement already satisfied: packaging in /usr/local/lib/python3.10/dist-packages (from geopandas) (23.2)
Requirement already satisfied: pandas>=1.1.0 in /usr/local/lib/python3.10/dist-packages (from geopandas) (1.5.3)
Requirement already satisfied: pyproj>=3.0.1 in /usr/local/lib/python3.10/dist-packages (from geopandas) (3.6.1)
Requirement already satisfied: shapely>=1.7.1 in /usr/local/lib/python3.10/dist-packages (from geopandas) (2.0.3)
Requirement already satisfied: attrs>=19.2.0 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (23.2.0)
Requirement already satisfied: certifi in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (2024.2.2)
Requirement already satisfied: click~=8.0 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (8.1.7)
Requirement already satisfied: click-plugins>=1.0 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (1.1.1)
Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (0.7.2)
Requirement already satisfied: six in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (1.16.0)
Requirement already satisfied: setuptools in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (67.7.2)
Requirement already satisfied: python-dateutil>=2.8.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.1.0->geopandas) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.1.0->geopandas) (2023.4)
Requirement already satisfied: numpy>=1.21.0 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.1.0->geopandas) (1.25.2)
Collecting rasterio
  Downloading rasterio-1.3.9-cp310-cp310-manylinux2014_x86_64.whl (20.6 MB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 20.6/20.6 MB 31.4 MB/s eta 0:00:00
Collecting affine (from rasterio)
  Downloading affine-2.4.0-py3-none-any.whl (15 kB)
Requirement already satisfied: attrs in /usr/local/lib/python3.10/dist-packages (from rasterio) (23.2.0)
Requirement already satisfied: certifi in /usr/local/lib/python3.10/dist-packages (from rasterio) (2024.2.2)
Requirement already satisfied: click>=4.0 in /usr/local/lib/python3.10/dist-packages (from rasterio) (8.1.7)
Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.10/dist-packages (from rasterio) (0.7.2)
Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (from rasterio) (1.25.2)
Collecting snuggs>=1.4.1 (from rasterio)
  Downloading snuggs-1.4.7-py3-none-any.whl (5.4 kB)
Requirement already satisfied: click-plugins in /usr/local/lib/python3.10/dist-packages (from rasterio) (1.1.1)
Requirement already satisfied: setuptools in /usr/local/lib/python3.10/dist-packages (from rasterio) (67.7.2)
Requirement already satisfied: pyparsing>=2.1.6 in /usr/local/lib/python3.10/dist-packages (from snuggs>=1.4.1->rasterio) (3.1.1)
Installing collected packages: snuggs, affine, rasterio
Successfully installed affine-2.4.0 rasterio-1.3.9 snuggs-1.4.7

Now, let's connect Drive to access the files that will be used:

In [ ]:
from google.colab import drive
drive.mount('/content/drive')

Mounted at /content/drive

We will import the libraries that we are going to use

In [ ]:
import geopandas as gpd
import numpy as np
import pandas as pd
import shapely
import math
import os
import rasterio
from rasterio.mask import mask
import json
from rasterio.plot import show
import matplotlib.pyplot as plt
import cv2
from skimage.morphology import skeletonize, convex_hull_image
from skimage.measure import label, regionprops
from shapely.geometry import shape
import shapely
from rasterio.features import shapes
from shapely.geometry import LineString

Then we set the path of the image we are going to use:

In [ ]:
path_img = '/content/drive/MyDrive/Datasets/Lines_planting/AOI_img.tif'

With rasterio we open the image and present it with matplotlib:

In [ ]:
src = rasterio.open(path_img)
In [ ]:
img = src.read()
In [ ]:
img = img.transpose([1,2,0])
In [ ]:
img = img.astype('uint8')
In [ ]:
plt.figure(figsize=[16,16])
plt.imshow(img)
Out[ ]:
<matplotlib.image.AxesImage at 0x7a4d6a4d8370>
No description has been provided for this image

In our RGB image we will separate the spectral bands and select the first spectral band.

In [ ]:
band = img[:,:,0]
In [ ]:
plt.figure(figsize=[16,16])
plt.imshow(band,cmap='RdYlGn')
plt.axis('off')
Out[ ]:
(-0.5, 7335.5, 7630.5, -0.5)
No description has been provided for this image

Let's look at the histogram distribution of this band:

In [ ]:
pixels = img[:,:,0].ravel()
pixels = pixels[pixels != 0]
In [ ]:
ax = plt.hist(pixels, bins = 256)
plt.show()
No description has been provided for this image
In [ ]:
band  = img[:,:,0].copy()
band = np.where(band == 0 , 255, band)

We select a value to divide the distribution into two classes:

In [ ]:
_,thresh1 = cv2.threshold(band,90,1,cv2.THRESH_BINARY_INV)
In [ ]:
plt.figure(figsize=[16,16])
plt.imshow(thresh1, cmap='Greys_r')
Out[ ]:
<matplotlib.image.AxesImage at 0x7a4d544daa70>
No description has been provided for this image

Now let's get the surrounding polygon of our region:

In [ ]:
hull = convex_hull_image(thresh1 == 1)
In [ ]:
plt.figure(figsize=[16,16])
plt.imshow(hull.astype('uint8'), cmap='Greys_r')
Out[ ]:
<matplotlib.image.AxesImage at 0x7a4d54369d20>
No description has been provided for this image

And convert the polygon image into a geodataframe that we will use later.

In [ ]:
shape_gen = ((shape(s), v) for s, v in shapes(hull.astype('uint8'), mask=hull, transform=src.transform))
Poly_gdf = gpd.GeoDataFrame(dict(zip(["geometry", "class"], zip(*shape_gen))), crs=src.crs)

Still with the image of the polygon, let's get the orientation of this polygon:

In [ ]:
for region in regionprops(hull.astype('uint8')):
  orientation = region.orientation
  angle_in_degrees = orientation * (180/np.pi)
  print(angle_in_degrees)

36.20655725077493

We will then apply a rotation to the RGB image using the angle we obtain:

In [ ]:
rows,cols,_ = img.shape
In [ ]:
M = cv2.getRotationMatrix2D(((cols-1)/2.0,(rows-1)/2.0),-32,1)
dst_img = cv2.warpAffine(img,M,(cols,rows))
In [ ]:
plt.figure(figsize=[16,16])
plt.imshow(dst_img)
plt.axis('off')
Out[ ]:
(-0.5, 7335.5, 7630.5, -0.5)
No description has been provided for this image

We will also apply rotation to the binary image:

In [ ]:
rows,cols, = thresh1.shape
In [ ]:
M = cv2.getRotationMatrix2D(((cols-1)/2.0,(rows-1)/2.0),-32,1)
dst = cv2.warpAffine(thresh1,M,(cols,rows))
In [ ]:
plt.figure(figsize=[16,16])
plt.imshow(dst)
plt.axis('off')
Out[ ]:
(-0.5, 7335.5, 7630.5, -0.5)
No description has been provided for this image

Let's now apply the morphological elements of opening and dilation using rectangular kernels to filter the lines in the image:

In [ ]:
kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (5, 25))
morph = cv2.morphologyEx(dst, cv2.MORPH_OPEN, kernel)

kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (50, 5))
morph = cv2.morphologyEx(morph, cv2.MORPH_DILATE, kernel)

# use morphology open to remove thin lines from dotted lines
kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (10,3))
morph = cv2.erode(morph, kernel,  iterations = 5)
In [ ]:
plt.figure(figsize=[16,16])
plt.imshow(morph)
plt.axis('off')
Out[ ]:
(-0.5, 7335.5, 7630.5, -0.5)
No description has been provided for this image

We apply skeletonization to obtain lines:

In [ ]:
from skimage.morphology import skeletonize, convex_hull_image
In [ ]:
sk = skeletonize(morph == 1)
In [ ]:
plt.figure(figsize=[16,16])
plt.imshow(sk,cmap='gray')
plt.axis('off')
Out[ ]:
(-0.5, 7335.5, 7630.5, -0.5)
No description has been provided for this image

Since these lines have a lot of vertical noise, we will filter the image with a sobel in the horizontal direction:

In [ ]:
sobely = cv2.Sobel(morph.astype('uint8'),cv2.CV_8U,0,1,ksize=5)
In [ ]:
plt.figure(figsize=[16,16])
plt.imshow(sobely,cmap='gray')
plt.axis('off')
Out[ ]:
(-0.5, 7335.5, 7630.5, -0.5)
No description has been provided for this image

With this we can return the image to its original position.

In [ ]:
rows,cols,_ = img.shape
In [ ]:
M = cv2.getRotationMatrix2D(((cols-1)/2.0,(rows-1)/2.0),32,1)
result_img = cv2.warpAffine(sobely,M,(cols,rows))
In [ ]:
plt.figure(figsize=[16,16])
plt.imshow(result_img)
plt.axis('off')
Out[ ]:
(-0.5, 7335.5, 7630.5, -0.5)
No description has been provided for this image

We can use Hough Lines to identify planting lines:

In [ ]:
from skimage.transform import hough_line, hough_line_peaks
from skimage.feature import canny
from skimage.draw import line
In [ ]:
tested_angles = np.linspace(-np.pi / 2, np.pi / 2, 360)

h, theta, d = hough_line(result_img, theta=tested_angles)
hpeaks = hough_line_peaks(h, theta, d, threshold=0.50 * h.max())
In [ ]:
for _, angle, dist in zip(*hpeaks):
  print(dist)

5305.0
7229.0
5743.0
4920.0
4791.0
8030.0
4640.0
4267.0
7743.0
3974.0
6269.0
6410.0
6123.0
6558.0
5590.0
5075.0
2518.0
3170.0
5890.0
7893.0
3846.0
7079.0
4128.0
5453.0
3472.0
7362.0
8175.0
3027.0
2661.0
4499.0
3699.0
3319.0
2090.0
6945.0
2370.0
6694.0
2896.0
2229.0
3485.0
3680.0
7524.0

Let's create a line geodataframe and plot it along with the RGB image:

In [ ]:
plt.figure(figsize=(16,16))
plt.imshow(img)
i = 0
crs_target = src.crs.to_dict()['init']
gdf_full = gpd.GeoDataFrame([])
lines = []
lines_id = []
row1, col1, _ = img.shape
for _, angle, dist in zip(*hpeaks):
    y0 = (dist - 0 * np.cos(angle)) / np.sin(angle)
    y1 = (dist - col1 * np.cos(angle)) / np.sin(angle)
    pt1 = (0,int(y0))
    pt2 = (col1,int(y1))
    retval, pt1_result, pt2_result = cv2.clipLine((0,0,col1,row1), pt1, pt2)
    plt.plot((pt1_result[0], pt2_result[0]), (pt1_result[1], pt2_result[1]), '-r')

    coordinates_pt1 = rasterio.transform.xy(src.transform, pt1_result[1], pt1_result[0])
    coordinates_pt2 = rasterio.transform.xy(src.transform, pt2_result[1], pt2_result[0])

    coords = [coordinates_pt1,coordinates_pt2]
    line = LineString(coords)
    lines.append(line)
    lines_id.append(i)
    i = i+1
gdf = gpd.GeoDataFrame(data=lines_id,geometry=lines)
gdf = gdf.rename(columns={0:"Line_Id"})
gdf = gdf.set_crs(crs_target)
gdf_full = pd.concat([gdf_full,gdf])
No description has been provided for this image

Looking at the lines geodataframe, we saw that the lines go outside the valid area of ​​the image:

In [ ]:
gdf_full.plot()
Out[ ]:
<Axes: >
No description has been provided for this image

Let's use the surrounding polygon geodataframe to clip out the parts of the lines that are not in the valid image area:

In [ ]:
Lines_cliped = gpd.overlay(gdf_full, Poly_gdf, how='intersection')
In [ ]:
Lines_cliped.plot()
Out[ ]:
<Axes: >
No description has been provided for this image

Finally we save it as a line shapefile:

In [ ]:
if not os.path.isdir('/content/shapefile'):
    os.mkdir('/content/shapefile')
In [ ]:
Lines_cliped.to_file("/content/shapefile/Lines.shp")
In [ ]:
!zip -r /content/shapefile.zip /content/shapefile
from google.colab import files
files.download("/content/shapefile.zip")

  adding: content/shapefile/ (stored 0%)
  adding: content/shapefile/Lines.prj (deflated 36%)
  adding: content/shapefile/Lines.shx (deflated 60%)
  adding: content/shapefile/Lines.cpg (stored 0%)
  adding: content/shapefile/Lines.shp (deflated 62%)
  adding: content/shapefile/Lines.dbf (deflated 91%)